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Abstract 

We investigate the stability of the square vortex lattice which has 
been recently observed in experiments on the borocarbide family of 
superconductors. Taking into account the tetragonal symmetry of 
these systems, we add fourfold symmetric fourth-derivative terms to 
the Ginzburg-Landau(GL) free energy. At H C 2 these terms may be 
treated perturbatively to lowest order to locate the transition from a 
distorted hexagonal to a square vortex lattice. We also solve for this 
phase boundary numerically in the strongly type-II limit, finding large 
corrections to the lowest-order perturbative results. We calculate the 
relative fourfold anisotropy for field in the xy plane to be 4.5% 
at the temperature, T C D , where the transition occurs at H C 2 for field 
along the z axis. This is to be compared to the 3.6% obtained in 
the perturbative calculation. Furthermore, we find that the phase 
boundary in the H — T phase diagram has positive slope near H C 2. 

1 Introduction 

The possibility of square vortex lattices in type-II superconductors was first 
proposed by Abrikosov M in his original paper on vortex lattices. However, it 
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was later shown that the hexagonal vortex lattice has a lower free energy for 
isotropic materials. Anisotropy in the effective mass or superfluid density 
tensors cause distortions of this hexagonal lattice simply in proportion to the 
anisotropy of the magnetic penetration length. Effects due to higher-order 
anisotropies (e.g., in cubic crystals) were seen in the 1970's for some low k 
superconductors such as PbTl, PbBi, and Nb. [[| In some cases, when the 
external magnetic field was applied along a fourfold symmetric axis of the 
crystal, the vortices formed a square lattice. When the field was applied 
along the threefold or twofold symmetric axis of the crystals, hexagonal or 
distorted hexagonal lattices were observed, respectively. In the case of Nb, for 
field along [001] the square vortex lattice transformed into a hexagonal lattice 
as the temperature or field increased. Takanaka [|J] introduced terms with 
cubic anisotropy into the phenomenological Ginzburg-Landau(GL) theory to 
explain the vortex lattice structures quite successfully. The extra terms arise 
from anisotropy of the Fermi surface. 

Recently, in ErNi 2 B 2 C, one of the borocarbide family of superconductors, 
Yaron, et al. |J and Eskildsen, et al. have used neutron scattering and Bit- 
ter decoration to observe the phase transition from a distorted hexagonal 
vortex lattice aligned with either the [100] or [010] directions to a square vor- 
tex lattice aligned with the [110] direction for increasing magnetic field along 
the [001] direction. This was the first observation of a square vortex lattice 
in a strongly type-II material with k» 1/V2. These experiments have mo- 
tivated further theoretical and experimental studies. Kogan et al. @ have 
analyzed the vortex lattice transition in borocarbides by taking into account 
the fourfold anisotropic corrections to the London model in the intermediate 
field region and T <C T c . This study revealed two kinds of phase transitions: 
at small fields (along the [001] direction), the distorted hexagonal lattice is 
aligned with a [110] direction. As the field increases, the lattice reorients 
along a [100] direction in a first-order phase transition. At a still higher 
field the square vortex lattice becomes stable. The latter phase transition is 
continuous. The square vortex lattice has also been observed in other boro- 
carbides such as LuNi 2 B 2 C and YNi 2 B 2 C (TIJ, using scanning tunnelling 
microscopy (STM) and neutron scattering, respectively. 

Vortex lattice structures in the high-temperature superconductors have 
also been studied. For low applied fields along the c axis (normal to the 
copper-oxide planes), the flux lattice can be seen using Bitter decoration. A 



distorted hexagonal lattice was seen for YBCO, [11] distorted as expected 
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from the in-plane (ab) effective mass anisotropy. But the STM pictures of 
the vortex lattice obtained at much higher field [12] are more consistent 
with a square lattice that is distorted by effective mass anisotropy. Thus it 
appears that the hexagonal to square lattice transition that we are studying 
may also occur in YBCO. Motivated by this, Affleck et al. |§ showed how 
fourfold anisotropic terms can arise in the London free energy due to a <i-wave 
symmetry of the superconducting order parameter. 

In this paper, we examine the structure of the vortex lattice near 
within a generalized Ginzburg-Landau (GL) free energy functional. We 
consider systems with tetragonal symmetry such as the borocarbides. The 
lowest-order term (ignoring effective mass anisotrpy) that reduces the sym- 
metry from isotropic down to tetragonal is fourth-order in the spatial deriva- 
tive and second-order in the order parameter. We discuss the three distinct 
terms that can be added to the GL free energy at this order. These new 
terms may be treated perturbatively near H c2 to first order, demonstrating 
that the square lattice occurs at high fields where the higher-order gradi- 
ent terms are more important, while the lattice distorts toward hexagonal 
at lower fields, in qualitative agreement with experiment. The phase transi- 
tion, which we treat only on a mean- field level, is an Ising-type critical point, 
with the square unit cell distorting into a rhombus in two possible ways by 
increasing the length of one of its diagonals and reducing the other. The per- 
turbation expansion has not been carried to higher order, so the quantitative 
accuracy of the first-order approximation had not previously been estimated. 
To investigate this we have numerically determined the phase diagram of this 
system, locating the upper critical field for various field orientations and the 
phase boundary separating the square and distorted hexagonal vortex lattice 
phases. We find that there are substantial (up to 50%) differences from the 
perturbative approximation, for example in the location of the phase bound- 
ary. Thus although the first-order perturbative approximation gives a good 
qualitative description of the phase diagram, quantitatively accurate results 
require a proper treatment of higher-order effects. 



2 Ginzburg-Landau model 

To study the lowest free energy pattern of the vortex lattice near H c2 for 
tetragonal crystals, we use the Ginzburg-Landau (GL) theory with added 
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fourfold symmetric fourth-order derivative terms. The full free energy func- 
tional that we minimize (ignoring thermal fluctuations) is 

F = F + F 4 . (1) 

F is the usual GL free energy: 

Fo = JdH*M 2 + ^ + ^m 2 + ^}. (2) 

For simplicity, we have taken an isotropic effective mass tensor. Tetragonal 
symmetry, of course, also allows uniaxial effective mass anisotropy, but this 
can be scaled out by rescaling the z axis and the xy components of the 
magnetic field. The only resulting difference is then in the h 2 field energy 
term in F , but this does not affect any of our results below. F 4 contains the 
additional terms allowed with tetragonal symmetry that are of fourth order 
in the x and y gauge-invariant derivatives and second order in the order 
parameter: 

+e a r(\(u 2 x - nj)^| 2 - KiLn, + ryi^l 2 ) 

c z 

The gauge-invariant derivative is 

h -> e* - 

IT = -V--A (4) 

I c 

We are interested in T < T c , so a is negative, r is defined to be 1 — 

J- c 

so it is positive. h(r) is the local magnetic field, h = V x A. contains 
two isotropic terms with coefficients e« and e&, and the anisotropic term with 
coefficient e a . The approximation of limiting the free energy to include only 
terms that are low-order in the gradient and the order parameter is valid for 
temperatures near T c . Only terms that are of fourth or higher order in the 
x and y gradients can reduce the rotational symmetry about the z-axis from 
circular down to tetragonal, which is why we add the fourth-order terms 
to the usual GL free energy. The anisotropic term may arise simply from 
fermi-surface anisotropy [|]] or from a d-wave nature of the pairing ||. 
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In the literature [§, 0], certain particular combinations of the isotropic 
and the anisotropic fourth order derivative terms in have been used. In 
our numerics (below), we will also use a specific convenient combination of 
them. For the theory as specified above to be stable, the fourth-derivative 
terms have to be non- negative, which requires > |e a |. In our notation, a 
term er|(II^ — n 2 )?/>| 2 corresponds to 

e e 3e 

€i, € a , £b — 2' 2' 2 ' 

A term er(|n^| 2 + |I1^| 2 ) || corresponds to 

ei,e a ,e b = T ,-;, T . (6) 

And the term 6x1(11^.11^ + n y IL E )'i/'| 2 that we use in our numerics corresponds 
to 

- _£ — (7\ 

In F4, the anisotropic term is the one that causes the vortex lattice to dis- 
tort and become square when the field is along the z axis and |e |r is large 
enough. The isotropic terms do not by themselves cause any distortion of 
the hexagonal vortex lattice, but they do change H c2 . 

For convenience, we introduce the following "dimensionless" units: the 
magnetic field is measured in units of H c2 = </> /(27r£ 2 ), the vector potential 



in units of H c2 £, } the length in units of the coherence length £ = h/ y2m*\a\, 
the order parameter in units of the zero-field order parameter magnitude: 



I ^00 1 — y M/A and the energy in units of H 2 ^ 3 /4w. Using these units, the 
free energy (D) is 



J rfV[-M 2 + + m 2 + ^ + e ,r|(n 2 + n 2 )^ 2 
+6 a r(|(n 2 , - n 2 )^| 2 - 1(11,11, + iwvf) 

+e b rh 2 M 2 ], (8) 



with 



n = — - a (9) 
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3 Perturbative results 



First we look at the vortex lattice behavior at H c2 by treating F± as a per- 
turbation to the usual GL theory, F , obtaining the free energy to first order 
in F4. Most of these results were first presented by Takanaka |4j], and then, 
much more recently, by De Wilde, et al. |J . 

Because of F 4 , the upper critical field depends on the field orientation. 
When the field is along the z axis, the upper critical field (in the absence of 
fluctuations) is given by 

# c2 [001] = ^(l-e t r-e b r + ...). (10) 

When the field is in the xy plane, 

H c2 [xyO] = ^(l--( ej r + ea rcos40) + ...), (11) 

where <fi is the angle between the field and the x axis. The relative xy-plane 
anisotropy of the upper critical field is thus 

AH c2 /H c2 = (F c2 [100]-# c2 [110])/F c2 [100] (12) 

= ~ + -' (13) 

Near H c2 the GL free energy density, including the first-order perturbative 
contributions from F4, may be simplified by the Abrikosov identities Jl], [|, 1^ . 
For the field along the z axis, this gives 

1 <h s > 2 2n 2 <h s > 2 



where ipL is the solution of the linearized GL equation at H c2 , including the 
terms from F±. h s is defined by h = H + h s , where H is the applied field; all 
fields are along the z axis. < ■ • • > means a spatial average over the vortex 
lattice, and B =< h > is the spatial average of the magnetic field. The free 
energy density (1|) can be written in terms of the Abrikosov ratio /3 a as 

, 2r R 2 (B ~ H c2 (e t T,e b T)) 2 

^ l + ^-lwW (15) 
_ 8/t 2 (e t r + t b r) - 
7 2k 2 -1 
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where 



Pa = 





| 4 > 


< 




2 >2 



Po - e a TReP 4 , (16) 



^ 2 " " <|Vol 2 > 2 2™ • (17) 

■00 is the solution to the linearized GL equations when = 0, and Il + is 
U x — illy. We obtain the stable pattern of the vortex lattice just below H C 2 
by minimizing the Abrikosov ratio, Pa- For regular vortex lattices with one 
vortex per unit cell, the lattice vectors are specified by (see Fig. 1) 

ai = ( T *-y), 

(The vortex lattice may also be rotated from this orientation by angle 7r/4; 
this is equivalent to changing the sign of e a .) p and /3 4 can be written 
explicitly using the reciprocal lattice vectors as 

Po = £exp[-^{(p 2 + g 2 )(l + & 2 )+2 M (l-& 2 )}], (18) 

p,q A0 
7T 2 

p,q 40 

exp [-|-{(p 2 + g 2 )(l + b 2 ) + 2pg(l - 6 2 )}], (19) 

where b = L x /L y . The sums run over all (positive, zero and negative) integer 
pairs. The square vortex lattice is b = 1 and the hexagonal lattice is b = 
or b = 1/ v 7 ^; other values of b we call distorted hexagonal lattices. 

For small positive |e |r the lowest free energy states for this system are 
distorted hexagonal lattices, and the square vortex lattices are local max- 
ima of the free energy. To calculate the stability of the square lattice the 
second derivatives at b = 1 are needed: they are d 2 Po/db 2 = —0.295 and 
d 2 RePi/ db 2 = 12.4. Thus, within the approximation of calculating Pa to 



7 



only first order in e a , at H c2 the appropriately-oriented square vortex lattice 
changes from a local maximum of the free energy to a global minimum at 



e a |r c u = 0.0238. 



(20) 



For larger values of |e |r the lowest free energy state of the vortex lattice 
is square. This critical value of |e a |r does not depend on k at first order in 
e a . This first-order approximation then says that the phase transition from 
distorted hexagonal to square vortex lattice at H c2 for the field along the z 
axis is at the temperature where the xy-plenae H c2 varies by about 3.6% from 
[100] to [110] field directions. 

One way to get a sense of how accurate this first-order in F4 approxima- 
tion is would be to calculate the next corrections (second-order in F4). This 
did not appear to be analytically tractable to us, so we have instead located 
the phase transition numerically, finding that there are indeed rather large 
corrections to this first-order approximation. 



At the order we are considering, in F4 there are three different terms added to 
the usual Ginzburg-Landau theory, and F contains k and H, which cannot 
be scaled away. In principle, one could examine the behavior in the full 
parameter space (H,K,€i,e a ,€b). However, we are interested in the phase 
transition to the square vortex lattice that is caused by the anisotropic term, 
e a . If we add only this term, it can be negative so the free energy then 
becomes unbounded below and the system is unstable at short wavelengths. 
To stabilize the system we set = \e a \. We also set e b to keep H c2 at its 
original value of (/> /(27r£ 2 ) for field along the z axis. For simplicity, we will 
only examine the large k limit, where the magnetic field is uniform. With 
all these constraints and simplifications, the system we studied numerically 
is, in dimensionless form: 



with er = 2ejT = — 2e a r > and for each value of er, e^r is chosen so that 
along the z axis H c2 = 1 in these units. In the first-order approximation 



4 Numerical Study 



F 




(21) 
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above, for field along the z axis the phase transition to the square vortex 
lattice at H C 2 occurs at er = 0.048. We find the higher-order corrections are 
substantial, and the actual transition in this system QZT| ) does not occur until 
eT c — 0.073. Presumably the precise result for e a r° will vary, depending on 
the values of k, and e& as well as other higher-order terms, all of which 
potentially differ between materials; we have not explored this dependence. 
Our basic message here is that the first-order in e a perturbative estimate 
is substantially off, so higher-order effects have to be treated properly to 
quantitatively estimate the phase diagram for any particular system. 

The stable vortex lattice for this system (^1|) with the field along the z axis 
is shown in Fig. 1. The area of the unit cell is dictated by the magnetic field 
because there is exactly one flux quantum per vortex. The only remaining 
freedom is the angle, 9, between the two nearest-neighbor vectors shown in 
Fig. 1. The square lattice is 9 = |, the hexagonal lattice is 9 = | (or 
other angles are the distorted hexagonal lattices. The minimum free energy 
vortex lattice for each 9 is an order parameter and current pattern that is 
symmetric under the operations that reverse the currents and then reflect in 
a plane that contains vortex lines and is parallel to the x or y axis. These 
symmetries mean we only have to determine the pattern over one quarter of 
the unit cell of the vortex lattice. We have chosen the quarter in the positive 
x and y quadrant (see Fig. 1). 

We discretize the system on a rectangular numerical grid aligned with the 
x and y axes. We use gauge-invariant variables only, namely the magnitude 

of the order parameter at each grid point, and the gauge- invariant phase 
differences between nearest neighbor grid points along the x direction, A x (p, 
and along the y direction, A y <p. These phase differences are not all indepen- 
dent, since their sum around an elementary plaquette of the grid must be 
equal to Hd x d y , where d x and d y are the grid spacings. We put one grid point 
at the precise core of the vortex at the origin in Fig. 1; at this point = 
so the phase is ill-defined and special care must be taken when treating the 
terms involving this point. 

All terms in the discretized GL free energy are symmetric so that the 
differences between the discretized expression and the continuum limit are 
even polynomials of d x and d y . We extrapolate to the continuum limit by 
reducing d x and d y together, keeping their ratio constant. Under these con- 
ditions the discretization errors are proportional to d x d y , the area of the grid 
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plaquette, for small d x and d y . For example, under our discretization 

.V, 



|(— -A x )^(x,y)\ 2 -> 
i 

^|[|^ + y,y)|exp(^) - |^-|,y)|exp(-^)]| 2 . (22) 
The boundary conditions on our quarter unit cell are 



[-0(0, 0)| 


= o, 




l#c,^)| 




— X 


= \«T 


A^z, 0) 


= o, 




A„0(O,j/) 


= o, 




A y 0(y,y) 


= o, 




A,0(x, ^) 


= -A x (f> 


— 

v 2 


A,0(x, ^) 


= -A y 





(23) 
(24) 

(25) 
(26) 

(27) 

*,^), (28) 
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x,^f); (29) 



we put grid points along all these boundaries. 

We begin with smooth initial conditions that satisfy all of the above 
constraints. Then we minimize the discretized version of the GL free energy 



I]) using the relaxation method [OJ] 

r e W = r l d _ T __ U (30) 

where T is a relaxation parameter which can be adjusted to avoid numerical 
instabilities and to hasten the convergence to the minimum. 

We have obtained the minimum of the discretized free energy for grid sizes 
(12 x 6), (16 x 8), (20 x 10), (24 x 12), (28 x 14), and(32 x 16) points for the 
square {6 = 90°) and two weakly distorted vortex lattices {6 = 88°, 86°). For 
each H and er the dependence of the free energy on 6 is extrapolated to the 
continuum limit. At H C 2 we also obtain the Abrikosov ratio /3a{^,6). The 
critical parameter value above which the minimum Abrikosov ratio is at the 
square lattice {6 = 90°) is er° = 0.0734. To determine the phase boundary 
below H c2 , we fix er and sweep the magnetic field to find the transition 
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Table V. Transition fields, H D , for a few er values. 



er 


H a /H c2 


0.0734 


1.000 


0.0800 


0.901 ± 0.002 


0.0860 


0.826 ± 0.001 


0.0930 


0.752 ± 0.001 


0.1000 


0.687 ± 0.001 



field, H n , above which the square vortex lattice is stable. Table |1] shows the 
transition fields for a few er values. Within our model, we certainly could 
follow the phase boundary to lower fields and temperatures, deeper into the 
vortex lattice phase. But the model, in the usual approach of Landau theory, 
has kept only certain terms that are low-order in an expansion in the order 
parameter, the gradients, and r. The neglected higher-order terms become 
more important as one moves away from T c to lower temperature and away 
from H c2 to lower fields. Thus the phase diagrams of real systems with 
different higher-order terms will increasingly diverge from the behavior we 
find in this specific model (|2T|) as one goes to lower temperature. This is why 
we did not bother to follow the phase boundary to lower temperature than 
shown in Table |I[ 

The physically meaningful phase diagram in the magnetic field, temper- 
ature plane is shown in Figure |2], with the field and temperature axes scaled 
by the values of H and r at the point r°) where the phase boundaries 
H D (T) and H c2 (T) meet. Although the point at which this happens is at 
roughly 50% larger er than the estimate from the perturbative results, our 
phase diagram when scaled this way is very close to the one obtained pertur- 
batively by De Wilde, et a/.,|| with the distorted hexagonal/square phase 
boundary meeting the upper critical field with a positive slope in the H, T 
plane (with this scaling, the slope of H n (T) is about 0.2, while the slope of 
H c2 (T) is -1.) 

Since, for field along the z axis, the critical er° at H c2 obtained numeri- 
cally is quite different from the perturbative estimate, we have also numeri- 
cally calculated for this system fl2"T|) the relative anisotropy of H c2 in the x-y 
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plane at er c n . We find that the dependence of H c2 on the angle between 
the field and the contains substantial harmonic content beyond the 

basic cos(40) term that is present in the first-order estimate, again indicating 
contributions beyond first-order in e. At er = er c D , the in-plane H c2 varia- 
tion from [100] to [110] directions is AH c2 /H c2 = 4.5%, which is again larger 
than the 3.6% estimated perturbatively. The in-plane anisotropy of H c2 (<j)) 
has recently been measured for LuNi 2 B 2 C, where it grows to near 10% at low 
temperature, indicating that the transition to the square lattice does meet 



the upper critical field for this material. 15 



5 Conclusion 

We have numerically determined the vortex lattice phase diagram for a gen- 
eralized Ginzburg-Landau model fl2~T|) that includes higher-order terms that 
reduce the rotational symmetry to tetragonal. Although the results are qual- 
itatively unchanged from the estimates obtained by treating the anisotropic 
term perturbatively to only first order, there are substantial quantitative 
differences. In particular, we find that for field along the z axis, the phase 
boundary in the field-temperature plane separating the distorted hexagonal 
and square vortex lattices meets the upper critical field line with a positive 
slope at a temperature T C D . If the field is instead put in the xy plane at this 
same temperature, the difference in H c2 between the [100] and [110] direc- 
tions is about 4.5%. Although this number must be a better estimate than 
the 3.6% obtained perturbatively, the full Ginzburg-Landau theory for any 
given material will contain other higher-order terms that may also have a 
small influence on this critical anisotropy value. Thus the precise value may 
vary somewhat away from our estimate. 
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Figure 1: Geometry of the vortex lattice. The unit cell is drawn by the 
dotted line. The black dots represent vortices. a 1: a 2 are the lattice vectors. 
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Figure 2: Phase diagram of our generalized Ginzburg-Landau model (|2l|) , 
scaled to the multicritical points (T C ,H = 0) and (T°,H°). The horizontal 
axis is the scaled T — T c : —t s = —t/t^, and the vertical axis is the scaled 
field: H s = H/H^. The open circles represent our data points fitted to a 
polynomial function (solid line). For comparison, De Wilde, et a/.'s result, 
scaled the same way, is drawn by the dashed line. 
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